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Abstract 

Genetic switching driven by noise is a fundamental cellular process in genetic regulatory networks. Quantitatively 
characterizing this switching and its fluctuation properties is a key problem in computational biology. With an 
autoregulatory dimer model as a specific example, we design a general methodology to quantitatively understand the 
metastability of gene regulatory system perturbed by intrinsic noise. Based on the large deviation theory, we develop new 
analytical techniques to describe and calculate the optimal transition paths between the on and off states. We also 
construct the global quasi-potential energy landscape for the dimer model. From the obtained quasi-potential, we can 
extract quantitative results such as the stationary distributions of mRNA, protein and dimer, the noise strength of the 
expression state, and the mean switching time starting from either stable state. In the final stage, we apply this procedure to 
a transcriptional cascades model. Our results suggest that the quasi-potential energy landscape and the proposed 
methodology are general to understand the metastability in other biological systems with intrinsic noise. 
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Introduction 

Stochasticity is an inherent property of living cells. Especially 
when the low copy number of species like the DNA and mRNA 
are taken into account, stochastic fluctuations can become 
significant and may qualitatively affect the behavior of the whole 
system [1,2]. To deal with these fluctuations, cells have evolved 
many mechanisms, of which genetic switch is a typical example. 
Cellular systems performing genetic switches usually consists of 
one positive feedback or double negative feedbacks [3,4]. 
Depending on the robustness of the feedbacks, cells can perform 
switches either spontaneously or on call [5]. 

Previous kinetic studies about cellular stochasticity have been 
formulated by using the generating function [6], system size 
expansion [7,8], large deviation theory (LDT) [9-13], or by 
employing WKB approximation to the chemical master equations 
(CMEs) [14,15], etc. However, only few of them take transcrip- 
tional noise into account explicitly. Some recent studies have 
shown that correlations between mRNA and protein levels do not 
always perform equally well in revealing genetic regulatory 
relationships [16,17], and the involvement of mRNA has a large 
effect on the switching times [18,19]. On the other hand, ever 
since Waddington's "epigenetic landscape" proposed in 1957 [20], 
the energy landscape have been widely used to provide intuitive 
illustration of the dynamics and evolution of genetic regulatory 
systems [1,11,21]. Thus it is important and desired to have an 
approach which can effectively determine the key features of a 
noisy gene expression system, such as constructing the corre- 



sponding "Waddington potential", identifying the transition paths 
between metastable states and computing the transition rates, etc. 

In this paper, we present a methodology to understand the 
metastability of the genetic switches in gene expression driven by the 
intrinsic noise based on LDT for Markov processes [22-24]. By 
explicitly taking mRNA noise into account, we obtain the most 
probable transition paths for off-to-on and on-to-off genetic switches 
through the geometric minimum action method (gMAM) [25]. 
Furthermore, we construct the global quasi-potential energy 
landscape, which is the rationalized version of the Waddington 
potential in this context. Based on the obtained quasi-potential, we 
obtain quantitative results for transition rates between metastable 
states and the intrinsic noise strength of gene expression state. We 
also consider the reduction of redundant dimensions if we are only 
interested in the energy landscape for partial components of the 
whole system. We successfully apply this methodology to a 
transcriptional cascades model. The relation between our and 
other approaches in literature is also discussed. From the authors' 
opinion, this framework is generally applicable for studying 
transitions between stable-saddle-stable fixed points with jump type 
noise generated by Gillespie type birth-death dynamics [26]. The 
quasi-potential energy landscape will be a powerful tool to unravel 
the metastable properties in more general biological processes. 

Models and Methods 

We illustrate our method through a specific two-state gene 
expression model in Fig. 1. We will refer to it as the "dimer 
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model" throughout the remainder of this paper. In this model, the 
gene at the active state transcribes mRNA with a much larger rate 
than it at the inactive state. Proteins translated from mRNA can 
aggregate into dimers that bind to the promotor site of the gene via 
a positive feedback. All of the processes are modeled as elementary 
reactions and all reaction rates are rescaled by protein decay rate 
(i.e. we will set y n = 1 unless stated otherwise). Here we assume that 
the reaction rates of dimer binding and dropping from DNA are 
much larger than the other reactions. 

Classical Methods and Issues 

The deterministic mean-field description of this dimer model 
through quasi-steady state approximation (QSSA) yields the 
ODEs. 



dM aoyo + akoD 
dt ~ y 0 + k 0 D 



„M, 



dN 

~di 



--bM- 



„N-2kiN 2 + 2 y] D, 



dD 
~dl 



--kyN 2 -y x D, 



where M, N and D are the mean number of the mRNA, protein 
and dimer respectively, and the parameters are shown in Fig. 1. 
This system has two stable fixed points and one saddle in 
physically reasonable parameter regime (see Text S 1 : 1) . These two 
stable fixed points correspond to the expressed and unexpressed 
states at which the copy number of proteins is at high or low state, 
respectively. With this deterministic description, once the system 
settles in one of its two attractive fixed points, it will stay there 
forever. However, in the presence of intrinsic noise, the system will 
fluctuate around its attractive fixed points and switch between 
these two metastable states on a large timescale. This type of 
switching among metastable states under small perturbations is 



To 



Figure 1. The autoregulatory dimer model with positive 
feedback. Promoter transitions are regulated by the dimerized 
transcription factor with rate ko and ro- a is the transcription rate of 
active promotor, with a very small transcription rate of inactive 
promotor ao- b is kinetic rate of translation, y m and y n are degradation 
rates of mRNA and protein, k\ and y x are the rates of dimerization and 
de-dimerization. All the processes are modeled as elementary reactions 
and all reaction rates are rescaled by the protein decay rate (i.e. y n = 1 
unless stated otherwise). 
doi:10.1371/journal.pone.0088167.g001 



exactly the rare events studied in the literature. The large 
deviation theory (LDT) is an appropriate tool to quantitatively 
describe the rare transitions [23,25,27]. Roughly it tells that when 
the system size K is sufficiently large, the probability that the 
trajectory of the stochastic dynamics X(i) stays in a small 8- 
neighborhood around a specific path (p can be given as 



sup 

0<t<T 



m 

K 



-<P(t) 



>S 



exp(-KI T (<p)), 



where It(<P) = Jo L(tp,ip)dt is called the rate functional. Thus the 
most probable transition path can be obtained by minimizing 
ItW) associated with the Lagrangian the L. Our task is to find the 
L for specific models. For Gillespie type birth-death dynamics, L 
has no closed form and only its dual Hamiltonian can be obtained 
in the large volume limit K— > oo, i.e. the number of all types of 
molecules goes to infinity. However, this approach encounters 
difficulty if we take the DNA switching into consideration since 
there is only one DNA copy in the considered model. Thus the 
straightforward utilization of the existed Hamiltonian in the large 
volume limit is invalid here. 

Large Deviation Theory 

To solve this issue, we develop the LDT direcdy for this specific 
system following the way in [28] with further extension. The 
biologically relevant choice of parameters in our model suggests 
the scaling }' 0 ,fl,ao ~ O(K), k\ ~ 0(1 /K) and the others are 0(1), 
where K is the system size which is usually chosen as the typical 
number of proteins in the expressed state. This is also the correct 
scaling under which the mean field limit of the CMEs gives the 
ODE system derived from QSSA (see Text SI:)). Define the 
rescaled concentration variable x = (x,y,z) where (x,y,z) = 
(m,n,d)/K and (m,n,d) is the state vector for the number of 
mRNAs, proteins and dimers. Correspondingly define 
y 0 = y 0 /K,a = a/K,ao=ao/K and k\=k\K to transform all the 
parameters to 0(1) magnitude. However for notational simplicity 
we will drop the tilde symbol on these parameters in the rest of the 
paper. It turns out that the Lagrangian of our model has the form. 



L(x,q)- 



inf (L\ (x,q,fi) - 



L 2 (x,fi)) 



(2) 



which combines the LDT result for large volume limit as shown in 
(1) and the Donsker-Varadhan type LDT result [22,29] for DNA 
fast switching. Here q = x resembles the velocity in classical 
mechanics, fl = (ji(0),{l(Y)) is a probabilistic 2-vector describes the 
residence distribution of DNA at the inactive or active state. 

Let us illustrate the construction of (2) via an intuitive way as 
follows. The net effect of the fast switching of DNA induces a 
residence distribution ft with components fi(0) and charac- 
terizing the probability of DNA staying at the inactive and active 
states, respectively. Whenever DNA is at the inactive or active 
state, we can apply the traditional LDT result in the large volume 
limit. This leads to the first part L\(x,q,fi) in (2). The second part 
describes how the visiting distribution induced by the random fast 
switching of DNA is close to a prescribed residence distribution fi 
given the current state x. This is exactiy what the Donsker- 
Varadhan type LDT gives [22,29], which is described by Li(x,ii). 
Since we are only interested in the LDT for the state variable x, 
the overall Lagrangian should be taken infimum with respect to all 
of the possible residence distributions fi. All of the statements will 
be made clear in the continued paragraphs. 
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Similar as the case in the large volume limit, it is not feasible to 
get the explicit form of the Lagrangian L{x,q) in general, but its 
dual Hamiltonian H{x,p) can be available, where p is the 
generalized momentum conjugate to x as in classical mechanics. 
They are connected via the Legendre transform. 

H(x,p)= SUp((q,p}-L(x,q)) 
1 



L 2 (x,f,)=-M £ ^^V), (8) 

" u(a) 

where u = (u(0),u(l)) is any 2-vector and C x is the infinitesimal 
generator for the DNA two-states jumping process at a given state 
x defined as 



supj <?,/»>- inf (Li(x,q,/i)+L 2 (x,/i)) 



'u(oy 




' -k 0 z 


k 0 z ' 




u(oy 






. Yo 


-Yo. 




Ml). 



The direct calculation shows that 



L 2 (x,fi) = ( sj fcozMO) - VToMl)) 



(9) 



(10) 



sup sup((q,py-L l (x,q,fi)-L 2 (x,fi)) 

1 M 



Combining Eqs. (3), (5), (6), (7) and (10), we obtain the final 
explicit Hamiltonian by optimization. 



sup suv((q,p}-Li(x,q,fi)-L 2 (x,/i)) 



= sup(Hi(x,p,/i)-L 2 (x,fi)). 

P 



(3) 



For the Gillespie type birth-death process with M reaction 
channels and the propensity function aAx) and stoichiometric 
vector Vj for j=\,2, . . . M, the existed LDT result gives the 
Hamiltonian [23]. 



H(x,p)=J2a J (x)( e v f-l). 



(4) 



Specifically in our model we have. 

Hi (x,p,fi)=ti(Q)H o[( (x,p) + n(\)H on {x,p) 
with the Hamiltonian 

H off (x,p) = a 0 (ei>x - 1) + A(x,p) 
when the DNA is at the inactive state and 

H on (x,p) = a(<? x - 1) + A(x,p) 



(5) 



(6) 



(7) 



when the DNA is at the active state. Here p = (p x ,Py,Pz) and 
A(x,p) = y m x(e-^ - \) + bx{eP>' - \) + y H y{ e -i>y - \) + k x y 1 (e- 1 ''y + ^ - 1) 
+ y\z(e 2p >'~ Pz — 1) corresponds to the part of the Hamiltonian for 
the gene expression, i.e. the transcription and translation 
processes. On the other hand, the famous Donsker-Varadhan 
LDT gives the Lagrangian 



H(x,p) = [as + a 0 (l -j)](e» - l)-(^/k 0 z(l-s) 



(11) 



where s= l/2 + ?/(2\A 2 +4) and t = (k 0 zy o y 1/2 [(a-a 0 )(eP x - 1) 
+ koZ — y 0 ]. This derivation can be easily extended to similar 
problems. 

It is worth noting that one can show the Hessian of H{x,p) with 
respect to p has the form. 



Vtr7 = 



d\ 


0 


0 


0 


d 2 + 4fi?3 — 2<5?3 


0 




d, 


dx 


0 0" 




0 


0 


d 2 0 


+ d] 


2 


0 


0 0 




-1 



(12) 



[0 2 -1], 



where d\ = 2(a - a 0 f / y' 'k 0 zy 0 (t 2 +A)eP x + (a 0 (l - s) + as)e Px + 
y m xe~ Px , d 2 = bxe p y +y n ye- p y and d 3 =k x y 2 e- 2p y +Pz + 
y l ze lp y~ Pz . Since d\,d 2 ,d^ are positive when x,y,z are positive, 
the decomposition in (12) shows that the Hessian is positive 
definite. This means the Hamiltonian is convex with respect to p in 
physically meaningful domain. It is quite different from that 
obtained by WKB asymptotics [15] (See Text SI: II). The 
convexity of the Hamiltonian is testified to be essential for the 
robustness and efficiency of the numerical algorithm [25] both 
theoretically and practically. It can be rigorously proved that it is a 
natural by-product from the LDT analysis. 

Quasi-Potential and Optimal Transition Path 

With the obtained LDT, we can get the optimal transition path 
through variational optimization. Moreover, we can define the 
local quasi-potential S(x\ Xq) with respect to a meta-stable state Xq 
as. 
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S(x:xo) = inf inf LUp,m)dt. 
r>o M <o)=x< ) , 9 <r)=iJo 



(13) 



S(x) = mm{S(x;x oa ) + AS,S(x;x 0 [f)}- 



(17) 



From classical mechanics, the local quasi-potential S(x; Xq) 
satisfies a steady-state Hamilton-Jacobi equation characterized by 
the Hamiltonian shown in Equation (11). 



H(x,V x S) = 0. 



(14) 



The LDT also confirms the equilibrium distribution of the 
system through a global quasi-potential function S(x). 



pWxexp{-B(x)}, 



(15) 



where the symbol x means the equality relation on a logarithmic 
scale, and S(x) can be obtained from its local version S(x; Xq) by a 
suitable sticking procedure which we will describe in the next 
subsection. This function S(x), which naturally serves as a 
rationalized version of the Waddington potential, is one main 
point of this article. 

The classical Hamilton-Jacobi theory enables one to solve the 
local quasi-potential S(x; Xq) satisfying (14) with variational 
methods. Here we employ the powerful geometric minimum 
action method (gMAM) proposed in [25] to compute S by 
minimizing the action functional with a prescribed Hamiltonian 
(11). The key idea of gMAM is essentially the Maupertuis principle 
in classical mechanics, which reformulates the action functional on 
the space of curves with intrinsic parameter, thus frees the time 
variable in the minimization process and still keeps its efficiency in 
high dimensions. This approach also resolves the issue of the 
singular boundary value problem by solving Hamilton's equations 
directly [15,30]. Specifically after each run of gMAM with fixed 
starting and ending points, one obtains the minimized action S 
and the corresponding optimal path. The readers may be referred 
to Text S1:IV for more details. 

It is worth asking whether the choice of the large parameter K 
affects the final results since any choice is artificial in practice. An 
affirmative answer is given in Text S 1:111 that only the scaling 
matters and the final systems are equivalent with respect to 
different choices of the large parameter K. 

Construction of Global Quasi-Potential Energy Landscape 

Based on the obtained local quasi-potential S(x) starting from 
the on and off states, we may construct the global quasi-potential 
energy landscape for genetic switching model by sticking them 
together. The system with only two metastable states and one 
saddle point, as our dimer model, can be handled conveniendy as 
the way shown below. The readers may refer to [24] for systematic 
methods of sticking the global quasi-potential for more complex 
systems. 

In our dimer model, we first compute the local quasi-potential 
S(x; Xo) starting from two metastable states Xrj = *on and Xrj = *off ■ 
We define S on = S(x s ; x on ) and S 0 ff = S(x s ; x 0 ff), where x s is the 
saddle point. Denote AS= \S on — S 0 s\. Suppose S on >S 0 (f, then 
the global quasi-potential S(x) is given by. 

S(x)= mm{S(x;x on ),S(x;x o[[ ) + AS}, (16) 
otherwise S(x) has the form 



It is not difficult to observe that sticking the two local quasi- 
potential via the linking saddle x s is the key point in this 
construction. 

In most cases, the considered system is in high dimensions while 
we are only interested in partial variables which is in low 
dimensions. This is also the case in our dimer model. Although the 
global quasi-potential S(x) is in three dimensions, we are mainly 
interested in its 2D configuration in the mRNA-Protein plane. So 
we need to reduce the redundant dimension z to obtain a 2-D 
potential S(x,y). We proceed with the following arguments. 

According to the LDT analysis (15), we obtain. 



P(x,y,z) x exp(-KS(x,y,z)). 
With the same reason, we have 

P(x,y)~exp(-KS(x,y)) 



(18) 



(19) 



where P{x,y) is the reduced distribution for mRNA and protein. 
By definition, this distribution is given by 

P(x, y ) = | P( x ,y,z)dz x | exp ( - KS(x,y,z))dz. (20) 

From the Laplace asymptotics [29] we get a simple reduction 
strategy. 



S(x,y) = min S(x,y,z). 



(21) 



This argument is general for any high dimensional situations 
and indeed it is also applied to the noise cascading model 
considered in our later text. 

Results 

Optimal Transition Path 

The large-deviation theory predicts that when events with little 
likelihood occurs, they will follow the optimal transition path 
which minimizes the action (13) with high probability. The 
probability of those paths deviated from the optimal one will decay 
exponentially. By choosing two stable states as the starting and 
ending points respectively, we obtained the switching path from 
either of the two states (see Fig. 2). For the convenience of 
visualization, we project the transition paths onto the mRNA- 
protein plane. 

Figure 2 shows clearly that when switch occurs, the trajectory 
prefers to be around the most probable path characterized by the 
Hamiltonian (11). The fact that the off-to-on and on-to-off paths 
are not identical agrees with the previous studies that the switching 
process is irreversible. The irreversibility is fundamental in 
chemical reaction kinetics due to the non-gradient nature of the 
considered system and can be considered as a form of hysteresis. 
However, in contrast with the previous study [11], our results 
indicate that when the noise level goes to zero, both optimal 
transition paths pass through the same bottleneck, i.e. the saddle 
point obtained from the corresponding deterministic model. This 
suggests that the saddle point has the lowest barrier height along its 
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number of mRNA number of mRNA number of mRNA 



Figure 2. Switching paths (A) from off to on state (purple solid curve) and (B) from on to off state (red solid curve) and MC 
simulations for both switching trajectories. We take the two stable fixed points in the deterministic dynamics as the starting and ending points. 
Darkness of the shading points represents the number of visits for reactive trajectories with smoothing. (C) Averaged switching trajectories from MC 
simulation. For each number of protein, we average in the mRNA dimension using probability as weight. Here the statistical results around each 
stable state is not shown because of the restrictions by our MC simulation algorithm (see Text SI:VI-A). The results are obtained from 1000 
independent long time MC simulations. The parameters here are ko = 1, y 0 = 50, ao = 0.4, a = 400, 6 = 40, y m = 10, y n = \,k\ =0.0002, and y 1 =2. 
doi:1 0.1 371 /journal.pone.00881 67.g002 



stable manifold and is in accordance with the energy landscape 
shown later. 

Global Quasi-potential Energy Landscape 

Applying the method of constructing global quasi-potential, we 
can compute the 2-D potential S(x,y) for our dimer model. The 
result is shown in Fig. 3. 

In Fig. 3, we observe that the on and off states correspond to 
two local minimum on the quasi-potential energy landscape, the 
saddle of the deterministic dynamical system exactly corresponds 
to the saddle point on the quasi-potential energy landscape too. 
The flatness along the mRNA direction keeps in good accordance 
with the large fluctuation observed in the reactive trajectories. 

To further characterize the switching path, we denote the first 
half (i.e. the part between the starting point and the saddle point) 
as the uphill path and the latter half as the downhill path. One 
may note that the transition path is also given by the Hamilton's 
equations x = V p H, p = — V X H. Therefore based on the fact 
H(x,O) = 0, we obtained V x H = 0 when p = 0. At the saddle point 
in any transition path, we have p = 0 [25] , and thus p = 0 along the 



A 



1200 



800 



400 




20 
15 

5 10 

1200 



whole downhill path. With this result we obtain the downhill 
equations x = V p H(x,0), which exactly corresponds to the 
corresponding deterministic dynamics. This fact explains that 
after climbing the saddle point the biological system relaxes to its 
attracting state fast without costing any action. This fact was also 
pointed out in [30-32]. 

On the other hand, the Hamilton-Jacobi theory also yields the 
uphill dynamics. 



x = V p H(x,V x S). 



(22) 



It is difficult to give a thorough understanding about the whole 
uphill path because of the general nonlinearity of H. However, an 
analysis around the critical points is instructive. Based on the fact 
p = 0 at critical points (i.e. the metastable states and saddle), we 
have by Taylor expansion. 



o o m 



n 

Figure 3. Quasipotential energy landscape of the whole genetic switching system with (A) two and (B) three dimensional view as 
well as switching paths between two stable fixed points. Each path passes through the saddle point. Here, the parameters are the same as in 
Fig. 2. 

doi:1 0.1 371 /journal.pone.00881 67.g003 
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xkV p H(x,0) + V 2 Mx,0)V x S. 



(23) 



JCASon 



(26) 



Recall that W p H(x,Q) corresponds to the deterministic mean 
field ODEs, the equation (23) is exactiy the uphill path of a 
chemical Langevin dynamics [33] (See more details in Text 
S1:VI). However, this chemical Langevin dynamics is not a 
straightforward generalization from the mean field ODEs like the 
usual large volume limit. Indeed, this reflects the specialty in our 
model setup which is related to the DNA fast switching. To see this 
more concretely, we have the approximated uphill dynamics for 
the x-component in our dimer model as. 



a 0 y 0 + ak 0 z 
k 0 z + y 0 

+ y m xS x . 



a 0 y 0 + ak 0 z 2{a- 



«o) 2 



k 0 z+y 0 k 0 z + y 0 



(24) 



The corresponding effective Langevin dynamics for the X- 
component reads. 



a 0 y 0 + ak 0 Z 
k 0 Z + y 0 



(25) 



where r\ x and rj 2 are two standard temporal Gaussian white noise 
with mean E)7,(f) = 0 and covariance i.(rj j (t)rjj(s)) = 8ijS(t — s), and 

a(z) = (a 0 y 0 + ak 0 z)/{k 0 z + y 0 ) + 2(a — a 0 ) 2 /(k 0 z + y 0 ). It is re- 
markable that besides the mean field transcription rate 
(aoyo + akoz)/(koz + y 0 ), we have an additional term 
2(a — do) 2 /(koz + y 0 ) in a. This is quite different from the usual 
chemical Langevin equations where for the j-th reaction the 
diffusion term sjaj(xj\j has the same form as its corresponding 
drift part aj(x)vj except a square root operation. The additional 
term in a makes that the fluctuation of the transcription is larger 
than that in gene expression processes without DNA switching, yet 
has the same mean field transcription rate 
(aoy 0 + akoz)/(koz + y 0 ). This observation coincides with the 
theoretical analysis in [34] although no feedback is considered 
there. 

The quasi-potential energy landscape not only provides the 
pictorial illustration for the dynamical transitions, it also contains 
many quantitative information to understand the metastability in 
genetic switching models. Once the global energy landscape is 
obtained, one can get the stationary distribution of the whole 
system via Eq. (15). Furthermore, it is very easy to calculate two 
main characteristic quantities used to describe a genetic switching 
system through only small amount of computational efforts. One is 
the transition rates, corresponding to the lowest barrier heights 
between two metastable states, and the other is the noise strength, 
corresponding to the steepness of quasi-potential around each 
metastable state. 

Global Property: Mean Switching Time (MST) 

The transition rate of switching systems is often characterized 
by Mean Switching Time (MST). We can compute the MST from 
either metastable state in the dimer model. For example, 
according to [24], the MST x on from on-to-off transition can be 
estimated from an asymptotic analysis. 



Here T on is a prefactor, and the quasi-potential energy barrier is 
AS on = So — S on , where So and S on are the action values at the 
saddle and on states, respectively. The result for the MST of off- to- 
on transition T 0 ff is similar. Although for one dimensional system 
the prefactor of MST can be obtained [31], there are no available 
results in high dimensions because of the geometry problem and 
the non-gradient nature of the system [35,36]. Fortunately, the 
prefactor varies slowly in many cases, therefore we can compare 
the MC simulations with the exponential time part and adjust the 
prefactor T on to fit the numerical results. 

The sensitivity of both MSTs, T on and T 0 ti, to the change of 
transcription rate a and mRNA decay rate y m are investigated and 
compared with MC simulations in Fig. 4A and Fig. 4B, where the 
prefactors of off-to-on and on-to-off transitions are estimated as 
2300 and 29, respectively. It can be observed that the MST is 
excellently predicted by Eq. (26) up to a slowly varying prefactor. 
And it is worth noting that when the MST becomes very large, the 
efficiency of the classical MC simulations gets extremely low while 
it is well kept in our approach. 

The positive feedbacks in genetic circuit usually provide cellular 
memory or all-or-none switch. The results in Fig. 4 reveal the 
robustness and sensitivity of the dimer model to the change of 
different kinetic parameters. In Fig. 4A and Fig. 4B, the promoter 
transition rates ko = l,yo = 50. When the transcription rate a 
increases from 360 to 420, the MST from on-to-off states T on 
increases exponentially, while the MST of off-to-on transition T 0 ff 
decreases slowly; it means the on-state becomes more stable while 
the off-state can still keep its stability. Therefore when the 
transcription rate a is increased, our genetic dimer circuit with 
positive feedback provides a stable cellular memory at the on state, 
but the off state remains its stability. Thus the system can not 
switch from the off state to the on state effectively. 

How to turn on the genetic switch? The results in Fig. 5 provide 
two possible effective choices. The first choice is to pose an 
additional source of the mRNA production, which we call as the 
trigger signal. We denote the additional mRNA production rate as 
Strigger- If St r igg er increases from 0 to 20, then MST of off-to-on 
transition T 0 ff will decrease exponentially and the genetic switch is 
turned on. We also show how the global energy landscape changes 
with different trigger rate Sagger = 5 in Fig. 5C and Strigger = 10 in 
Fig. 5D. It is evident to see from the figures that when the trigger 
signal increases, the barrier height from off-state to on-state 
decreases, which is in accordance with the MST of off-to-on 
transition t d iy. The relevant biological switch can be found in the 
start point of budding yeast cell cycle process, where the additional 
trigger signal in Gl cyclin Cln3 causes the activation of Gl 
transcription factor SBF and MBF [37]. The second possible 
choice is to decrease the degradation rate of protein y„. When the 
degradation rate of protein y n decreases from 1 to 0.5 in Fig. 5B, 
MST from off-to-on transition T c ff will also decrease exponentially 
from 10 6 to 10 4 and turn on the genetic switch. This is the case in 
Bacillus subtilis, where the transitions into competent state is caused 
by decreasing the degradation rate of protein ComK [38]. All the 
other relative parameters are listed in the caption of Figures. 

Furthermore, we calculate the MST of both from off to on and 
on to off states in the different promoter transition rates. We show 
the results with the fast rates (ko = 5, y 0 = 250) in Fig. 4C and slow 
rates (ko = 0.2, yo = 10) in Fig. 4D. These results indicate that the 
system with slow promoter transition rates tend to have short MST 
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Figure 4. The mean switching time (MST) and quasipotential energy landscape as a function of parameters. (A) and (B): MST as a 
function of transcription rate a. Promoter transition rates kg = 1, y 0 = 50, the gMAM results with numerical prefactor of off-to-on transition (red solid 
line) and on-to-off transition (blue dashed line), compared with MC simulations (•) and (o), respectively. (C) and (D): The gMAM results with different 
promoter transition rates of off-to-on transition (red) and on-to-off transition (blue), where solid line with &o = l, 7o = 50 is same as (A) and (B), the 
faster transition rate in dashed line with &o = 5,y 0 = 250, the slower transition rate in dotted line with ko=0.2, y 0 = 10. Other parameters are 
ao = 0.4, 6 = 40, y„ = 1, ki =0.0002, y 1 =2; in (A,C), y m = 10 and (B,D) n = 400. 
doi:10.1371/journal.pone.0088167.g004 



of off-to-on transition (red dashed lines) while fast rates lead to long 
MST of off-to-on transition (red dash-dotted lines). This is due to 
the reason that faster promoter transition rates lead to smaller 
mRNA and protein noise strength (see Fig. 6 for more detailed 
information). Ignoring the difference of the mechanism of initial 
transcription between prokaryotes and eukaryotes, in the simple 
case, the faster promoter transition rates correspond to the gene 
expression process in prokaryotes, and the slower promoter 
transition rates correspond to the slow chromatin remodeling 
process in eukaryotic case [39]. The results suggest that 
prokaryotes may have stronger cellular memory than eukaryotes. 

Local Property: Fluctuation Around Stable States 

Another quantitative information that quasi-potential energy 
landscape can provide is the noise strength of stable states. Here 
we use the coefficient of variation (CV, i.e. the standard deviation 
over the mean) to measure the strength of fluctuation instead of 
the Fano factor, for the system here has positive feedback thus 
deviates far from Poisson statistics. Notice that the stationary 
distribution p(x) x exp{ — KS(x)}, we can expand S(x) in the 
vicinity of high stable state x on up to second order thus get the 
Gaussian approximation. 



1 1 _, T 

P(x)- 3 T exp{--(x-fi)I, (x-fi) }. (21) 

(27t)2|2|2 1 

Here, x = (x,y,z),fi = (x 0 „,y 0 „,z 0 „), S = (S,y) 3x3 , and |Z| is the 
determinant of matrix Z. Eq. (27) holds only in the vicinity of the 
on state with standard deviations a,„ = {K(S" V S", — S"?)/|X|)*, 

a„ = (K(S^ : - S'£)/\L\f and a d = (K(S^ y - S%)/\L\f. 
With the ff,„, a „ and Gd above, we can easily obtain the CV as 
shown in Fig. 6. 

Figure 6 demonstrates that when the average expression levels 
increase, the noise strength of mRNA and protein decreases in our 
positive feedback model. The fluctuation of mRNA is usually 
larger than that of protein. Furthermore, the noise level with slow 
promoter transition rates is almost always larger than the one with 
fast promoter transition rates. This is in accordance with the 
results of MST that the system with long MST has small noise and 
vise versa. The inconsistent portion between analytical and 
simulation results (the left part of the line with slow promoter 
transition rates in Fig. 6B) is due to the inapplicability of Eq. (27) 
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Figure 5. The mean switching time (MST) of off-to-on transition as a function of (A) trigger signal strength that transcribes mRNA at 



constant rate and (B) degradation rate of protein y„. (C) and (D): Quasipotential energy landscape with different trigger strength. S b 



= 5 in 



(C), and S u 



10 in (D). Other parameters are a 0 = 0.4, a = 400, 6 = 40, y, = 1, k x =0.0002, 



10; 



■■ 1 in (A,C,D), and S tri 



= 0 in (B). 



doi:1 0.1 371 /journal.pone.00881 67.g005 

during the low barrier crossing process for the on state. More 
details may be referred to Text S1:V. 

Application in Transcriptional Cascades 

To further illustrate the power of quasi-potential energy 
landscape and the abundant quantitative information it contains, 
we apply our methodology to a transcriptional cascades model 
based on the previous work of S. Hooshangi et al. [40]. In their 
work, S. Hooshangi et al. synthesized transcriptional cascades 
comprised of one, two, and three repression layers and analyzed 
the sensitivity and noise propagation as a function of network 
complexity. They used different concentrations of anhydrotetra- 
cycline (aTc) as inducer and measured the fluorescence intensities 
of protein eyfp (the last layer of each cascade) by the flow 
cytometer. 

Here we simplify the 3-layer cascades as x— *y\ H )>2 H J'3, 
where x denotes the concentration of aTc as inducing signal and 
j/l , _V2, J'3 denote the output of proteins in different layers 
respectively. Then we directly construct the quasi-potential energy 
landscape for each layer and obtained the normalized probability 
distribution of the output to certain signal x from Eq. (15). The 
dose response curves to increasing signal x axe shown in Fig. 7, 
which are consistent well with the previous experimental results. 
Further more, two features of transcriptional cascades can be 
observed. Firstly, the more layers the transcriptional cascades 
have, the sharper the response curves are (as the Hill coefficient of 
the 3-layer cascades is 2.00, 3.15 and 4.08 respectively). Thus the 
sensitivity is increased in the cascades. Secondly, the fluctuation of 
output can be described by the spreading width of its distribution, 
so more layers of cascades amplify the cell-cell variability (see Fig. 



S2). In short, when a cascade has more layers, its response curve 
gets steeper with a wider probability distribution and thus larger 
fluctuations. The straightforward calculation of CV based on Fig. 7 
has been done and it agrees well with the MC simulations (see 
Text S1:VI and Figure S2). 

Limitations of The Study, Open Questions, and Future 
Works 

We have already illustrated a general methodology based on 
LDT to quantitatively understand the metastability in gene 
expression processes perturbed by the intrinsic noise and applied 
it to a dimer auto-regulatory circuit model. It is clear that this 
methodology can be extended to more general systems, provide 
one can explicrdy write down the Hamiltonian of the system. If all 
of the considered species have relatively large numbers, the 
Hamiltonian is simply the Eq. (4). For the case where the large 
volume limit fails to be true, our method is also applicable under 
an additional assumption that the low copy number of species 
reach their stationary distribution much faster than the others. 
This is the situation that we treat DNA in our dimer model. 
However, we would like to mention the limitations of our work, 
which of course motivates us for future studies. 

The main limitations or the corresponding open questions can 
be summarized into the following three aspects: 

1 . The case where the large volume limit and the fast switching 
mechanism are both invalid. This prevents us to construct the 
LDT for the considered system. Thus there is no Hamiltonian 
and the current methodology fails. How to quantitatively study 
such systems and define the proper Waddington energy 
landscape is an issue. 
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Figure 6. The coefficient of variation (CV) versus mean number of (A) mRNA and (B) protein induced by varying transcription rate a 
with different promoter transition rates. The lines and discrete dots correspond to analytical results and MC simulations, respectively. The 
results with fast promoter transition rates are shown in blue dash-dotted line and □, medium rates in red solid line and o, and slow rates in magenta 
dashed line and 0- The parameters here are the same as in Fig. 4 (A,C). 



2. The curse of dimensionality. When the problem is extended to 
high dimensions, although the computation of optimal 
transition path and the analysis of MST and CV can be 
achieved with a reasonable cost, constructing the whole global 
energy landscape is not feasible in general. However as we have 
argued before, the whole energy landscape is even not 
necessary since we are only interested in its configuration for 
partial components. How to extract these information directly 
from the Hamiltonian by smart utilization of Eq. (21) is what 
we are trying to do. 

3. Transition rate formula in the high dimensional case. Despite 
the transition rate formula, i.e. the Arrhenius type formula, for 
the equilibrium models are well developed [41], there is no 
complete answer for the non-equilibrium case. The rate with 
form R x exp( — KAS) has long been established in [24], but 
the prefactor is not known. In one dimensional case, partial 
result is given [31]. But its high dimensional form is still an 
open question. 

To understand the transition behavior for more general 
biological systems driven by noise, the above open problems 
should be overcome in the future studies. 



Conclusion and Discussion 

In this paper, we have presented a methodology to construct 
the quasi-potential energy landscape of genetic switching system 
while explicitly taking mRNA noise into account. This global 
potential, which is a rationalized version of Waddington 
potential, can provide a quantitative tool to understand the 
metastability in more general biological processes with intrinsic 
noise. The results also provide some insights in gene-expression 
switching circuit with positive feedback, especially the robustness 
and sensitivity of the genetic switching system under different 
promotor transition rates. 

For the connection with previous general methodology in 
literature [11], we focus more on the energy landscape and 
metastability properties for systems with intrinsic noise. Although 
one can principally compute the stationary distribution by solving 
a steady state chemical master equation on a truncated domain, 
our approach sufficiently utilizes the special structure of the 
system. Indeed, the global quasi-potential S(x) employed in this 
paper connects with the potential defined in [1 1] through 
S(x)= — limA-^oo K~ l \np ss (x), which is independent of K. 
The reason we can do this is simply because the system size K 
is large enough here. 




10' 10° 10' 10 2 10"' 10° 10' 10* 10"' 10° 10' 10 ! 



Figure 7. The dose response curves and probability distribution of the output protein in the 3-layer cascades (denoted by r, .VT.ri) 
as a function of inducing signal x, 1 -layer in (A), 2-layer (B) and 3-layer (C). The probability distribution can be directly obtained from Eq. 
(15) after normalization. The Hill coefficient for each cascade is fitted as 2.00, 3.15 and 4.08 respectively. 
doi:10.1371/journal.pone.0088167.g007 
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With regard to the WKB and reduction approach in [15], we 
obtain a convex Hamiltonian based on rigorous mathematical 
analysis and explicitly take mRNA noise into account through the 
gMAM method. The convexity proves to be essential for the 
computational efficiency and robustness. Our derivations can be 
also easily extended to similar problems. 

Overall, the quasi-potential energy landscape and the proposed 
methodology can serve as a useful tool to explore the gene 
expression process with intrinsic noise. Further developments such 
as high dimensionality issue and its applications to other biological 
systems like complex cellular decision making and the develop- 
mental process of cells are deserved to be investigated. The 
biological meaning of optimal transition path and transition states 
remains to be uncovered in the future studies. 

Supporting Information 

Figure SI The network design of three synthetic 
transcriptional cascades. 

(EPS) 
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